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ABSTRACT 

Using high-resolution 3-D and 2-D (axisymmetric) hydrodynamic simulations in spherical geometry, 
we study the evolution of cool cluster cores heated by feedback-driven bipolar active galactic nuclei 
(AGN) jets. Condensation of cold gas, and the consequent enhanced accretion, is required for AGN 
feedback to balance radiative cooling with reasonable efficiencies, and to match the observed cool 
core properties. A feedback efficiency (mechanical luminosity « eM acc c 2 ; where M acc is the mass 
accretion rate at 1 kpc) as small as 6 x 10 -5 is sufficient to reduce the cooling/accretion rate by ~ 10 
compared to a pure cooling flow in clusters (with M 2 oo < 7 x 10 14 M 0 ). This value is much smaller 
compared to the ones considered earlier, and is consistent with the jet efficiency and the fact that only 
a small fraction of gas at 1 kpc is accreted on to the supermassive black hole (SMBH). The feedback 
efficiency in earlier works was so high that the cluster core reached equilibrium in a hot state without 
much precipitation, unlike what is observed in cool-core clusters. We find hysteresis cycles in all our 
simulations with cold mode feedback: condensation of cold gas when the ratio of the cooling-time 
to the free-fall time {t coo \/ts) is < 10 leads to a sudden enhancement in the accretion rate; a large 
accretion rate causes strong jets and overheating of the hot ICM such that Cooi/iff > 10; further 
condensation of cold gas is suppressed and the accretion rate falls, leading to slow cooling of the core 
and condensation of cold gas, restarting the cycle. Therefore, there is a spread in core properties, 
such as the jet power, accretion rate, for the same value of core entropy or t coo i/ts. A fewer number 
of cycles are observed for higher efficiencies and for lower mass halos because the core is overheated 
to a longer cooling time. The 3-D simulations show the formation of a few-kpc scale, rotationally- 
supported, massive (~ 10 11 M 0 ) cold gas torus. Since the torus gas is not accreted on to the SMBH, 
it is largely decoupled from the feedback cycle. The radially dominant cold gas (T < 5 x 10 4 K; 
v r > |^|) consists of fast cold gas uplifted by AGN jets and freely-infalling cold gas condensing out 
of the core. The radially dominant cold gas extends out to 25 kpc for the fiducial run (halo mass 
7 x 10 14 M 0 and feedback efficiency 6 x 10 -5 ), with the average mass inflow rate dominating the 
outflow rate by a factor of « 2. We compare our simulation results with recent observations. 

Subject headings: galaxies: clusters: intracluster medium - galaxies: halos - galaxies: jets 


1. INTRODUCTION 

Majority of baryons in galaxy clusters are in the form 
of a hot plasma known as the intracluster medium (ICM). 
In absence of cooling and heating, the ICM is expected to 
follow self-similar profiles for density, temperature, etc., 
irrespective of the halo mass (Kaiser 1986, 1991; see 
also the review by Voit 2005). However, self-similarity 
is not observed in either groups or clusters (e.g., Ponman 
et al. 1999; Balogh, Babul, & Patton 1999; Babul et 
al. 2002). Moreover, the core cooling times in about a 
third of clusters is shorter than 1 Gyr, much shorter than 
their age (~Hubble time; e.g., Cavagnolo et al. 2009; 
Pratt et al. 2009). Thus, we expect cooling to shape 
the distribution of baryons in these cool-core clusters. 

The existence of cool cores with short cooling times in 
a good fraction of galaxy clusters is a long-standing puz¬ 
zle. According to the classical cooling flow model, cluster 
cores with such short cooling times were expected to cool 
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catastrophically and to fuel star formation at a rate of 
100 — 1000 M 0 yr _1 (e.g., Fabian 1994; Lewis et al. 
2000). However, cooling, dropout, and star formation 
at these high rates are never seen in cluster cores (e.g., 
Edge 2001; Peterson et al. 2003; O’Dea et al. 2008). 
This means that some source(s) of heating is(are) able 
to replenish the core cooling losses, thereby preventing 
runaway cooling and star-formation. 

While there are potential heat sources, such as the 
kinetic energy of in-falling galaxies and sub-halos (e.g., 
Dekel & Birnboim 2008), thermal conduction from the 
hotter outskirts (e.g., Voigt & Fabian 2004; Voit 2011), 
a globally stable mechanism, which increases rapidly with 
an increasing hot gas density in the core, is required 
to prevent catastrophic cooling. Observations of several 
cool-core clusters by Chandra and XMM-Newton have 
uncovered AGN-jet-driven X-ray cavities, whose mechan¬ 
ical power is enough to balance radiative cooling in the 
core (e.g., Bohringer et al. 2002; Birzan et al. 2004; 
McNamara & Nulsen 2007). The AGN jets are powered 
by the accretion of the cooling ICM on to the supermas- 
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sive black hole (SMBH) at the center of the dominant 
cluster galaxy. Thus, more cooling/accretion leads to an 
enhanced jet power and ICM heating, closing a feedback 
loop that prevents runaway cooling in the core. 

AGN feedback has been long-suspected to play a role 
in self-regulating the ICM (e.g., Binney & Tabor 1995; 
Ciotti & Ostriker 2001; Soker et al. 2001; Babul et al. 

2002; McCarthy et al. 2008), but a clear picture has 
emerged only recently. While AGN feedback should pro¬ 
vide feedback heating in cluster cores (as it is enhanced 
with ICM cooling), it is not obvious if, for reasonable 
parameters, AGN heating can keep pace with cooling 
that increases rapidly with an increasing core density. 
Moreover, the dense core gas is expected to be highly 
susceptible to fragmentation, leading to the formation 
of a multiphase medium consisting of cold dense clouds 
condensing from the hot diffuse ICM itself. Pizzolato & 
Soker (2005) suggest that AGN outbursts that result in 
the heating of the cluster cores are due to the infall and 
accretion of these cold clumps. 

The importance of cold gas precipitation/feedback has 
also been been highlighted by several recent observations. 
The fact that there is some multiphase-cooling/star- 
formation, albeit at a much smaller rate than predicted 
by the cooling flow estimate (Soker et al. 2001), ties well 
with the idea of a small fraction of the thermally un¬ 
stable core gas cooling to the stable atomic and molec¬ 
ular temperatures. A lot of this cold gas is expected to 
form stars, but some should be accreted on to the cen¬ 
tral SMBH. Reservoirs of atomic (e.g., Crawford et al. 
1999; McDonald et al. 2011a; Werner et al. 2014) and 
molecular gas (e.g., Donahue et al. 2000; Edge 2001; Sa¬ 
lome et al. 2006; Russell et al. 2014; O’Sullivan et al. 

2015), both extended and centrally concentrated, and 
ongoing star formation (e.g., Bildfell et al. 2008; Hicks, 
Mushotzky, & Donahue 2010; McDonald et al. 2011b) 
are observed in a lot of cool-core clusters. Additionally, 
powerful radio jets/bubbles observed in most cool-core 
clusters (Cavagnolo et al. 2008; Mittal et al. 2009) can 
be interpreted as a signature of kinetic feedback due to 
cold gas accretion on to the SMBH. 

Since cool cores are in rough global thermal balance 
(i.e., the cooling rate minus the heating rate is smaller 
than just the radiative cooling rate), the existence of cold 
gas in cluster cores can be understood as a consequence 
of local thermal instability in a weakly stratified atmo¬ 
sphere (McCourt et al. 2012; Sharma et al. 2012a; Singh 
& Sharma 2015). The idealized simulations, which im¬ 
pose global thermal equilibrium in the ICM, show that 
the nonlinear evolution of local thermal instability leads 
to in-situ condensation of cold gas only if the ratio of 
the cooling time and the free-fall time (Aooi/te) is < 10 
(Sharma et al. 2012a). 

This model has the attractive feature that once the 
local thermal instability sets in and the cold gas begins 
to condenses out of the dense ICM, it typically falls freely 
toward the center. Some of the infalling cold gas has 
sufficiently low angular momentum to be accreted by the 
SMBH, resulting in the cold phase mass accretion rate 
onto SMBHs that can exceed the hot/Bondi accretion 
rate by a factor ~ 10 — 100 (Gaspari et al. 2013; Sharma 
et al. 2012a). This enhanced accretion rate in the cold 
phase can explain both the global thermal balance in 
cluster cores and the general lack of massive cooling flows 


in almost all cool-core clusters whereas, the hot-mode 
(Bondi) accretion rate appears inadequate by orders of 
magnitude (e.g., McNamara, Rohanizadegan, & Nulsen 
2011 ). 

In detail, the precipitation of the cold gas, followed by 
a sudden increase in the accretion rate onto the SMBHs, 
leads to an increase in jet/cavity power and (slight) over¬ 
heating of the core. The core expands and as the ratio 
t c ooi/tff rises above the threshold value of t coo \/tg=l0 , 
the gas is no longer prone to condensation. The accretion 
rate drops, as does the jet power. The core cools slowly 
and the whole cycle starts again when t coo i/tg< 10. The 
frequency of heating/cooling cycles depend on jet effi¬ 
ciency and the halo mass. These features of the cold 
feedback model are verified in our numerical simulations. 

In fact, the simple criterion of t coo i/tg< 10 for the on¬ 
set of local thermal instability is expected to be generic - 
applicable not only to the intracluster medium (ICM) but 
also the intragroup medium (IGrM) and the circumgalac- 
tic medium (CGM) of all galaxies, including the Milky 
way (Sharma et al. 2012b; Voit et al. 2015b). This, in 
turn, has far-reaching implications for providing a com¬ 
mon framework for understanding the the breaking of 
self-similarity in the properties of hot gas across the hi¬ 
erarchy, from galaxies to groups to clusters, the presence 
of multi-phase gas in group and clusters cores, and the 
detection of cold gas in galaxies at distances of ~ 100 kpc 
(e.g., Werk et al. 2014). In fact, recent more realistic 
AGN jet feedback simulations show that cold gas con¬ 
densation begins when the t coo \/tg< 10 condition is met, 
and two distinct cold gas structures emerge: extended 
cold filaments which go out 10s of kpc; and a few-kpc 
rotationally-supported cold torus (Gaspari et al. 2012; 
Li & Bryan 2014a,b). This dichotomy in cold gas dis¬ 
tribution is also seen in observations (e.g., McDonald et 
al. 2011a). 

Now that the theoretical models are satisfactorily able 
to describe the basic state of the ICM in cool cluster 
cores, and since observations of cold gas and jets/cavities 
are rapidly accumulating, it is ripe to make detailed 
comparisons between observations and numerical simula¬ 
tions. We also aim to investigate the similarities and dif¬ 
ferences in cold gas and jet/bubble properties as a func¬ 
tion of the halo mass and feedback efficiency. 

In this paper, we focus on cool-core clusters and have 
carried out 3-D and 2-D (axisymmetric) simulations of 
the interaction of feedback-driven AGN jets with the 
ICM over cosmological timescales, varying the halo mass 
and the feedback efficiency. The 3-D simulations, which 
should correspond more closely to reality, show the for¬ 
mation of a cold, massive, angular-momentum-supported 
torus, as seen in previous works (Gaspari et al. 2012; Li 
& Bryan 2014a,b). This massive cold torus is decoupled 
from the AGN feedback cycle, which is governed by the 
low angular momentum, radially-dominant (|iy| > \v$\, 
Vr/vtf, is the radial/azimuthal component of the velocity) 
in-falling cold gas. Angular-momentum-supported gas is 
absent in 2-D simulations because of axisymmetry and 
the absence of rotation in the initial state (stochastic 
angular momentum can be generated in 3-D because of 
d/dcj) terms in the angular momentum equation). How¬ 
ever, 2-D simulations are useful for two reasons: first, 
they show similar behavior to 3-D simulations, if we only 
consider the radially-dominant (|ty| > lu^j) cold gas; sec- 
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ond, they are much cheaper to run for long timescales, 
and thus are useful to do parameter scans in halo mass 
and accretion efficiency. 

Compared to previous works (Gaspari et al. 2012; Li 
& Bryan 2014a,b), we have carried out simulations with 
smaller feedback jet efficiencies. We find that a feedback 
efficiency as low as 6 x 10 -5 (ratio of the input jet power 
and M acc c 2 , where M acc is the accretion rate measured at 
1 kpc) is sufficient to reduce the mass accretion/cooling 
rate by a factor of about 10 compared to the cooling flow 
value in groups and clusters. Such a low feedback effi¬ 
ciency fits in nicely with the observations which suggest 
that only a small fraction (~ 0.01) of the available gas is 
accreted by the SMBH (e.g., Loewenstein et al. 2001), 
and with the estimate of jet efficiency (~ 0.001 — 0.01) 
with respect to the SMBH accretion rate (e.g., Benson 
& Babul 2009). Moreover, jets in SMBHs are observed 
predominantly when the accretion rate is < 0.01 the Ed¬ 
dington value (e.g., Narayan & Yi 1995; Merloni et al. 
2003); i.e., < 0.22 M 0 yr” 1 for a 10 9 M 0 SMBH. The 
expected mass accretion rate on to the SMBH in our 
simulations (^0.01 times M acc hi Table 1) satisfies this 
constraint. 

We have analyzed the velocity-radius distribution of 
the cold gas in our simulations to compare with re¬ 
cent ALMA and Herschel observations of cold gas struc¬ 
ture and kinematics in galaxy/cluster cores (e.g., Mc¬ 
Namara et al. 2014; David et al. 2014; Werner et 
al. 2014). Our simulations help in interpreting obser¬ 
vations of cold gas outflows and inflows at scales > 10 
kpc, and the rotationally-supported cold torus at scales 
< 5 kpc. In our simulations, the fast (> 500 km s -1 ) 
atomic/molecular outflows are uplifted by the outgoing 
AGN jet. The slower (< 300 km s -1 ) infall of cold gas is 
due to condensation in the dense core. The cold gas in 
the rotationally supported torus is at the local circular 
velocity (~ 200 km s -1 ). 

Our paper is organized as follows. In section 2 we 
present the numerical setup, in particular our implemen¬ 
tation of mass and kinetic energy injection due to AGN 
jets. Section 3 presents the key results from our 3-D and 
2-D simulations, a comparison of 3-D vs. 2-D, and the 
impact of parameters such as feedback efficiency and halo 
mass on our results. In section 4 we discuss our results 
and compare with previous simulations and observations, 
and we conclude with a brief summary in section 5. 

2. NUMERICAL SETUP & GOVERNING EQUATIONS 

We modify the ZEUS-MP code, a widely-used finite- 
difference MHD code (Hayes et al. 2006), to simulate 
cooling and AGN feedback cycles in galaxy clusters. We 
solve the standard hydrodynamic equations using spher¬ 
ical (r, 9, <fi) coordinates, with cooling, external gravity, 
and mass and momentum source terms due to AGN feed¬ 
back: 


% + v-(ptO = s p , 

dpv „ , 

-Qj- +V.(pv®v) = 


(1) 

-Vp - pV$ + SpVjetr, (2) 

( 3 ) 


where p is the mass density, v is the fluid velocity, 
p = (7 — l)e is the pressure (e is the internal energy 
density and 7 = 5/3 is the adiabatic index), A(T) is 
the temperature-dependent cooling function, n e (n*) is 
the electron (ion) number density given by p/[p e (i) m p] 
(p e = 1.18 and pi = 1.3 are the mean molecular weights 
per electron and per ion, respectively, for the ICM with 
a third solar metallicity). For the cooling function, we 
use a fit proposed in Sharma, Parrish, & Quataert 2010 
(their Eq. 12 and solid line in their Fig. 1) with a stable 
phase at 10 4 K. 

In addition to the terms shown in Eqs. 1-3, the code 
uses the standard explicit artificial viscosity, and has 
implicit diffusion associated with the numerical scheme 
(Stone & Norman 1992). In addition to the standard 
non-linear viscosity, we use the linear viscosity, as rec¬ 
ommended by Hayes et al. (2006) for strong shocks (see 
their Appendix B3.2). 

We use a fixed external NFW gravitational potential 
$(r) due to the dark matter halo (Navarro et al. 1996); 

Tf 1— GM 2 qo l n (l + C200 r / r 200) /.■, 

[r> ~ r [ln(l + C 2 oo)-C2oo/(l + C2oo)]’ U 

where M 2 oo (^200) is the characteristic halo mass (radius) 
and C200 = f 2 oo/r s is the concentration parameter; the 
dark matter density within r 2 00 is 200 times the critical 
density of the universe and r s is the scale radius. In this 
paper we focus on cluster and massive cluster runs with 
M200 = 7 x 10 14 M 0 and 1.8 x 10 15 M 0 , respectively, 
and adopt c 2 oo = 4.7 for all models. 

We include the source terms S p for mass and S p Vj e tr 
for the radial momentum to drive AGN jets (uj et is the 
velocity which the jet matter is put in). 1 These source 
terms and the cooling term (in Eq. 3) are applied in an 
operator-split fashion. The mass and momentum source 
terms are approximated forward in time and centered in 
space. The cooling term is applied using a semi-implicit 
method described in Eq. 7 of McCourt et al. (2012). 

Our simulations do not include physical processes like 
star formation and supernova feedback. Star formation 
may deplete some of the cold gas available in the cores 
(see Li et al. 2015), but this is unlikely to change our 
results for a realistic model of star formation. Super¬ 
nova feedback is energetically subdominant compared to 
AGN feedback, and cannot realistically suppress cluster 
cooling flows (e.g., Saro et al. 2006). We only include 
the most relevant physical processes, namely cooling and 
AGN jet feedback, in our present simulations. 

2.1. Jet Implementation 

Jets are implemented in the active domain by adding 
mass and momentum source terms as shown in Eqs. 1 
& 2. The source terms are negligible outside a small 
biconical region centered at the origin around 9 = 0, n, 
mimicking mass and momentum injection by fast bipolar 
AGN jets. 

The density source term is implemented as 
S p {r,9) = NMj et i/j{r,9), 

1 We have also carried out narrow-jet simulations with momen¬ 
tum injection in the vertical [£] direction, but do not find much 
difference from our runs with momentum injection in the radial 
[r; see Eq. 2] direction. 


n e riiA(T), 
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where Mj et is the single-jet mass loading rate, 


ip(r, 0) 


2 + tanh 



+ tanh 


x 


1 + tanh 




(9jet + 0 - 7T 


(Je 


(5) 


describes the spatial distribution of the source term 
which falls smoothly to zero outside the small biconical 
jet region of radius rj et and half-opening angle 0j e t- We 
smooth the jet source terms in space because the Kelvin- 
Helmholtz instability is known to be suppressed due to 
numerical diffusion in a fast flow if the shear layer is unre¬ 
solved (e.g., Robertson et al. 2010). The normalization 
factor 


Znrfcti 1 - cosset) 


ensures that the total mass added due to jets per unit 
time is 2Mj e t- All our simulations use the following jet 
parameters: ay = 0.05 kpc, (9j et = 7 t/6 , and ay = 0.05. 
The jet source region with an opening angle of 30 degrees 
may sound large but we get similar results with narrower 
jets. Also, the fast jet extends well beyond the source 
region and is much narrower (c.f. third panel in Figure 
1). The jet radius rj e t is scaled with the halo mass; i.e., 


rj e t = 2 kpc 


/ M 2 oq \ 1 ^ 3 

\7 x 10 14 M 0 y 


The jet mass-loading rate is calculated from the current 
mass accretion rate (M acc ) evaluated at the inner radial 
boundary such that the increase in the jet kinetic energy 
is a fixed fraction of the energy released via accretion; 
i.e., 

Afjet'ufet = eMaccC 2 . (6) 

We choose the jet velocity i>j e t = 3 x 10 4 km s _1 (0.1c; c 
is the speed of light); such fast velocities are seen in X- 
ray observations of small-scale outflows in radio galaxies 
(Tombesi et al. 2010). The jet efficiency (e; our fiducial 
value is 6 x 10 -5 ) accounts for both the fraction of the 
in-falling mass at the inner boundary (at 1 kpc for the 
cluster runs) that is accreted by the SMBH and for the 
fraction of accretion energy that is channeled into the jet 
kinetic energy. Our results are insensitive to a reasonable 
variation in jet parameters (uj e t, 7>t, 0j e t, ay, oy), but 
depend on the jet efficiency (e). 

Like Gaspari et al. (2012), the jet energy is injected 
only in the form of kinetic energy; we do not add a ther¬ 
mal energy source term corresponding to the jet. We note 
that Li & Bryan (2014b) have shown that the core evolu¬ 
tion does not depend sensitively on the manner in which 
the feedback energy is partitioned into kinetic or ther¬ 
mal form. Another difference from previous approaches, 
which use few grid points to inject jet mass/energy, is 
that our jet injection region is well-resolved. 


2.2. Grid, Initial & boundary conditions 

Most AGN feedback simulations evolved for cosmolog¬ 
ical timescales (e.g., Gaspari et al. 2012; Li & Bryan 
2014a) use Cartesian grids with mesh refinement. How¬ 
ever, we use spherical coordinates with a logarithmically 


spaced grid in radius, and equal spacing in 8 and <f>. The 
advantage of a spherical coordinate system is that it gives 
fine resolution at smaller scales without a complex al¬ 
gorithm. Perhaps more importantly, a spherical set up 
allows for 2-D axisymmetric simulations which are much 
faster and capture a lot (but not all) of essential physics. 

We perform our simulations in spherical coordinates 
with O<0<7r, 0 < <f> < 27r, and r m j n < r < r max , with 

7’[min,max] = [1, 200] k P C x 1 O i”m 0 ) ‘ 

According to self similar scaling, we have scaled all length 
scales in our simulations (inner/outer radii 7- m i n /r max , 

1 /3 

t- 200, jet radius r jet ) as M 2( j 0 . 

We apply outflow boundary conditions (gas is allowed 
to leave the computational domain but prevented from 
entering it) at the inner radial boundary. We fix the 
density and pressure at the outer radial boundary to the 
initial value and prevent gas from leaving or entering 
through the outer boundary. Reflective boundary condi¬ 
tions are applied in 9 (with the sign of v^ flipped) and 
periodic boundary conditions are used in cj). We noticed 
that cold gas has a tendency to artificially ‘stick’ at the 
8 boundaries (mainly in 2-D axisymmetric simulations) 
for our reflective boundary conditions. This cold gas can 
lead to an unphysically large accretion rate close to the 
poles, and hence artificially enhanced feedback heating 
(Eq. 6). Therefore, we exclude 8 grid-points at each 
pole when calculating the mass accretion rate; these ex¬ 
cluded angles correspond to only 0.5% of the total solid 
angle for 128 grid points in the 8 direction. All our diag¬ 
nostics (M acc , entropy profiles, etc.) also exclude these 
small solid angles close to the poles. 

The resolution for 3 — D runs is 256 x 128 x 32 and 
for 2 — D runs is 512 x 256. Since we use a logarithmic 
grid in the radial direction, the resolution for 256 (512) 
grid points in the radial direction corresponds to a good 
resolution of A r/r = 0.02 (0.01). The minimum resolu¬ 
tion in the radial direction for the fiducial 3D (2D) run is 
ss 0.02 (0.01) kpc. For such a resolution our integrated 
quantities (mass accretion rate, jet power, cold gas mass, 
etc.) are converged. 

We focus on simulations of a galaxy cluster with 
M 2 qo = 7 x 10 14 M 0 but with different parameters 
such as feedback efficiency. For comparison we also car¬ 
ried out simulations for a massive cluster with M 200 = 
1.8 x 10 15 Mq. The initial conditions are the same as in 
Sharma et al. (2012a); i.e., we assume the initial entropy 

profile (K = Tkev/77-e^ 3 ; 7k e v is the ICM temperature in 
keV and n e is the electron number density) of the form 

K(r) = K 0 + K lm (^y\ (7) 

as suggested by Cavagnolo et al. (2009). 2 For our cluster 
runs, we set Kq = 10 keV cm 2 and A'100 = 110 keV cm 2 
at the start (as in Sharma et al. 2012a). We assume self¬ 
similar behavior scaling with M200 (Kaiser 1986) to set 

2 Whether an entropy core exists is debated (Panagoulia, Fabian, 
& Sanders 2014), but our results are insensitive to our initial 
conditions. Our ICM profiles change with time and reach a quasi¬ 
steady state which may or may not have an entropy core. 
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the initial entropy profile for our massive cluster runs (i.e. 
we assume K 0 = 19 keV cm 2 and K W q = 210 keV cm 2 ; 
c.f. McCarthy et al. 2008). Except for early transients, 
our results are independent of the precise choice of the 
initial values of Kq and Aboo- 

The outer electron number density is fixed to be n e = 
0.0015 cm -3 . Given the entropy profile and the density 
at the outer radius, we can solve for the hydrostatic den¬ 
sity and pressure profiles in an NFW potential (Eq. 4). 
We introduce small (maximum overdensity is 0.3) iso- 
baric density perturbations on top of the smooth density 
(for details, see Sharma et al. 2012a). 

3. RESULTS 

In this section we describe the key results from our 
simulations. Table 1 lists our runs. We begin with the 
results from our fiducial 3-D cluster run (C6m5D3 in Ta¬ 
ble 1). We show that the 1-D profiles of density, entropy, 
etc. are consistent with observations. We highlight the 
cycles of cooling and AGN jet feedback, and the spa¬ 
tial and velocity distribution of the cold gas. We show 
that there are three components in cold (T < 5 x 10 4 
K) gas distribution: a massive, centrally-concentrated, 
rotationally-supported torus; spatially extended and fast 
(> 500 km s" 1 ) outflows correlated with jets; and slower 
(< 300 km s" 1 ) in-falling cold gas that condenses out 
because of local thermal instability. Then we compare 
the results from our 3-D and 2-D axisymmetric simula¬ 
tions. We also explore the dependence of our results on 
the halo mass and the jet efficiency. 

3.1. The fiducial 3-D run 

We experimented with different values of jet efficiencies 
(e; Eq. 6) in our 3-D cluster ( M 200 = 7 x 10 14 M 0 ) 
simulations, and found that the average mass accretion 
rate for e = 6 x 10 -5 was about 10% of a pure cooling flow 
(see Table 1). Therefore, we choose this as our fiducial 
value, which is smaller compared to the values chosen 
by some recent works (Gaspari et al. 2012; Li & Bryan 
2014a,b), but is consistent with observational constraints 
(e.g., O’Dea et al. 2008). Our fiducial value should be 
considered as the smallest efficiency that is required to 
prevent a cooling flow in a cluster (this critical efficiency 
depends on the halo mass, as we shall see later). 

The minimum ratio of the cooling time (t coo i = 
3nfcsT/[2n e riiA]) and the local free-fall time (tg = 
[2 r/g} 1 / 2 = [2r 3 /GM(< r)] 1 / 2 ) is 7 for the initial ICM; 
this ratio (Aooi/tff) is a good diagnostic of the state of the 
cluster core in rough thermal balance. Since the initial 
condition is in hydrostatic equilibrium, there is negligi¬ 
ble accretion through the inner boundary, and therefore 
there is no jet injection. However, after a cooling time 
in the core (rs 200 Myr) there is a rise in the accretion 
rate across the inner boundary (M acc ), and hence in jet 
momentum injection (Eq. 6). The jet powers a bubble 
that heats the core and raises Aooi/tff, keeping the mass 
accretion rate well below the cooling flow value (c.f. top 
panel of Fig. 9). After this time the cluster core is in a 
state of average global thermal balance between radiative 
cooling and feedback heating via AGN jets. 

3.1.1. Jets, bubbles & multiphase gas 


Figures 1 show the snapshots (r — 6 plane at (f> = 0) of 
pressure, density, and temperature at different times for 
our fiducial 3-D run. The X-ray emitting ICM plasma is 
quite distinct from the dense cold (10 4 K) gas and from 
the low-density jet/bubble. The cold gas accreting on 
to SMBH gives rise to AGN jets. Before a cooling time 
(0.2 Gyr) there are no signs of cooling and jets. After a 
cooling time, accretion rate through the inner boundary 
(at 1 kpc) increases and bipolar jets are launched (0.29 
Gyr). The jets are not perfectly symmetric, as they are 
shaped by the presence of cold gas in their way. The 
inhomogeneities in the ICM enhances mixing with (and 
stirring of) the ICM core, resulting in effectiveness of our 
jets even with a low efficiency. 

Jets are fast in the injection region but become slow, 
buoyant, and almost in pressure balance with the ICM 
(compare the upper and middle panels of Fig. 1 ) because 
of turbulent drag and sweeping up of the ICM. In absence 
of further power injection, the bubbles are detached from 
the jets and rise buoyantly and mix with the ICM at 10s 
of kpc scales (3.15 Gyr in Fig. 1). Most of the cold gas 
is very centrally concentrated (within 10 kpc), but does 
condense out at larger radii, although never beyond 30 
kpc. 

As jets plough through the dense cold gas clouds, for¬ 
ward shock moves ahead of these clouds after partially 
disrupting them. The collision results in a reverse shock 
and a huge back-flow of hot jet material which mixes with 
the cooler ICM, driving the core entropy to higher val¬ 
ues. These back-flows and mixing are mainly responsible 
for heating the cluster core. 

3.1.2. Radial profiles 

Before discussing the detailed kinematics of cold gas 
and jet cycles, we show in Figure 2 the 1-D profiles of im¬ 
portant thermodynamic quantities (entropy [Tkev/^e^ 3 ] 1 
tcooi/tfi, Tkev) as a function of radius for the fiducial 
3-D run. In addition to the instantaneous profiles (at 1, 
2, 3, 4 Gyr), the median profile and spread about it are 
shown. The median is calculated for the entropy mea¬ 
sured at 20 kpc (roughly the core size) and all the profiles 
with entropy within one standard deviation at the same 
radius are shown in grey. 

The spread in quantities outside ~ 20 kpc is quite 
small, but increases toward the center because mul¬ 
tiphase cooling (leading to density spikes) and strong 
jet feedback (leading to overheating) are most effective 
within the core. The density at 1 Gyr is peaked toward 
the center, indicating that the cluster core is in a cooling 
phase. The spikes in density at 3 Gyr have corresponding 
spikes in entropy and t coo \/ts profiles, but not as promi¬ 
nent in the temperature profile. The temperature fluc¬ 
tuations are rather modest compared to fluctuations in 
other quantities because of dropout and adiabatic cool¬ 
ing. Temperature profiles show a general increase with 
radius, as seen in observations. 

There is a large spread in entropy toward lower values 
about the median at radii < 10 kpc (top-left panel in 
Figure 2). This is because there are short-lived cooling 
events during which the entropy in the core decreases sig¬ 
nificantly (simultaneously, density increases and t coo i/tg 
decreases). On the other hand, the increase in the core 
entropy is smaller but lasts for a cooling time, which is 
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Fig. 1.— The pressure (upper panel), electron number density (middle panel), and temperature (lower panel) contour plots (R — z plane 
at cj) = 0) in the core at different times for the 3-D fiducial run. The density is cut-off at the maximum and the minimum contour level 
shown. The low-density bubbles/cavities are not symmetric and there are signatures of mixing in the core. The left panel corresponds to 
a time just before a cooling time in the core. The second panel from left shows cold gas dredged up by the outgoing jets. The rightmost 
panel shows in-falling extended cold clouds. The pressure maps show the weak outer shock, but the bubbles/cavities so prominent in 
the density/temperature plot are indiscernible in the pressure map, implying that the bubbles are in pressure equilibrium and buoyant. 
Also notice the outward-propagating sound waves in the two middle pressure panels in which the jet is active. The in-falling/rotationally- 
supported cold gas has a much lower temperature and pressure than the hot phase. The arrows in the temperature plots denote the 
projected gas velocity unit vectors. 
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TABLE 1 
Table of runs 


Label 

dim. 

N r X Ng X 

min. resolution 
(kpc) 

M200 

(M 0 ) 

jet efficiency (e) 

<M acc >l 

(Meyr- 1 ) 

(-^facc,hot)/ (^facc,cold) 

jet duty 
eyelett (%) 

C6m5D3^ 

3 

256 x 128 x 32 

0.02 

7 X 10 14 

6 

X 10 -5 

25.1 (244.2)1 

0.2 (0.06) 

59.6 

C5m4D3 

3 

256 x 128 x 32 

0.02 

7 x 10 14 

5 

x 10~ 4 

7 

0.78 

82 

Clm2-D3 

3 

256 x 128 x 32 

0.02 

7 x 10 14 


0.01 

1.9 

1.3 

99.8 

C6m5D2t 

2 

512 x 256 x 1 

0.01 

7 x 10 14 

6 

x 10~ 5 

23.5 (170)L 

0.23 (0.19) 

64.8 

C6m6D2 

2 

512 x 256 x 1 

0.01 

7 x 10 14 

6 

CO 

1 

O 

X 

153 

0.27 

47.3 

Clm2-D2 

2 

512 x 256 x 1 

0.01 

7 x 10 14 


0.01 

0.77 

14.8 

99.8 

Clm4D2 

2 

512 x 256 X 1 

0.01 

7 x 10 14 


10- 4 

13.7 

0.3 

63.8 

C5m4D2 

2 

512 x 256 x 1 

0.01 

7 x 10 14 

5 

x 10~ 4 

5.1 

1.6 

72.9 

M6m6D2 

2 

512 x 256 x 1 

0.014 

1.8 X 10 15 

6 

CO 

1 

0 

X 

293 (299) 

0.3 (0.5) 

0.0 

M6m5D2 

2 

512 x 256 x 1 

0.014 

1.8 x 10 15 

6 

x 10~ 5 

77.7 

0.58 

50.1 

Mlm4D2 

2 

512 x 256 x 1 

0.014 

1.8 x 10 15 


10~ 4 

48.18 

0.62 

47.1 

M5m4D2 

2 

512 x 256 x 1 

0.014 

1.8 x 10 15 

5 

x 10~ 4 

18.7 

2.9 

63.8 

Mlm2-D2 

2 

512 x 256 x 1 

0.014 

1.8 x 10 15 


0.01 

8.1 

3.7 x 10 5 

99.8 


Notes 

‘C’ in the label stands for a cluster (M200 = 7 x 10 14 Mq) and ‘M’ for a massive cluster (M200 = 1.8 x 10 15 Mq). Label C6m5D3 indicates 
that it is a cluster run in 3-D with an efficiency e = 6x 10 —5 (Eq. 6). 
iThe fiducial 3-D and 2-D runs. 

1 Angular brackets denote time average over the full run. The quantities in brackets denote values for a pure cooling flow (M c f). Note that 
8 grid points close to the poles are excluded when calculating the accretion rates. 

Jet duty cycle is defined as the fraction of total time for which the jet power is > 10 40 erg s —1 . 
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Fig. 2.— X-ray emissivity-weighted (considering only 0.5-8 keV gas) 1-D profiles of important thermodynamic variables as a function of 
radius. Snapshots at 1, 2, 3, 4 Gyr are shown. Various quantities are obtained by combining 1-D profiles of density and pressure. The 

median and standard deviation (cr) of entropy (K = Ti ce v/ n e //3 ) at 20 kpc are calculated. Various profiles corresponding to the median 
entropy at 20 kpc (14 keV cm 2 ) are shown in different panels (black lines with ‘+’). Thick grey lines show the profiles for which the entropy 
at 20 kpc is within 1 — cr of its median value. 
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longer in this state. This behavior is generic, fairly insen¬ 
sitive to parameters such as the feedback efficiency and 
the halo mass. 

3.1.3. The cold torus 

While Figure 1 shows that cold gas can be dredged up 
by AGN jets (second panel; see also Revaz, Combes, & 
Salome 2008; Pope et al. 2010) and can also condense 
out of the ICM at large scales (fourth panel), majority 
of cold gas is at very small scales (< 5 kpc) in the form 
of an angular-momentum supported cold torus. Figure 3 
shows the zoomed-in density snapshots in the equatorial 
(0 = 7t/2) plane at different times; the arrows show the 
projection of velocity unit vectors. As the cluster evolves 
the cold gas, condensing out of the hot ICM, gains an¬ 
gular momentum from jet-driven turbulence. Because 
of a significant angular velocity, an angular momentum 
barrier forms and cold gas circularizes at small radii. 

Unlike Li & Bryan (2014b), our cold torus is dynamic 
in nature as AGN jets disrupt it time and again, but it 
reforms due to cooling. Figure 3 shows the evolution of 
the torus at various stages of the simulation. The top- 
left panel of Figure 3 shows the cluster center at 0.5 Gyr. 
Small cold gas clouds are accumulating in the core after 
the first active AGN phase. At 1.3 Gyr, cold gas accret¬ 
ing through the inner boundary has an anti-clockwise ro¬ 
tational sense. At 1.98 Gyr, cold gas (and the hot gas out 
of which it condenses) is rotating clockwise. Jet activity 
leading up to this phase has reversed the azimuthal ve¬ 
locity of the cold gas. At all times after this the dynamic 
cold gas torus rotates in a clockwise sense, essentially be¬ 
cause the mass (and angular momentum) in the rotating 
torus is much larger than the newly condensing cold gas. 

The torus gets disrupted due to jet activity but forms 
again quickly. The snapshots at 2.4 and 2.45 Gyr show 
that the inner region is covered by the very hot/dilute 
jet material. If the jets were rapidly changing direction 
as argued by Babul et al. (2013), we would in fact 
expect the cold gas torus to be occasionally disrupted 
by the jets. In the present simulations, however, this 
behavior is an artifact of our feedback prescription; we 
scale the jet power with the instantaneous mass inflow 
rate through the inner boundary (see Eq. 6). Even small 
oscillations of the cold torus can sometimes lead to a large 
instantaneous mass inflow through the inner boundary 
and hence an explosive jet event. The reassuring fact is 
that these explosive ‘events’ are rare and the jet material 
is quickly mixed with the ICM after these. In reality, 
most of the cold gas in the torus will be consumed by 
star-formation. Only the low angular momentum cold 
gas that circularizes closer in (< 100 pc) can be accreted 
by the SMBH at a short enough timescale. 

A cold torus forms in all our 3D cluster simulations 
with different efficiencies. However extended cold gas is 
lacking at late times in simulations with high jet efficien¬ 
cies. Li & Bryan (2014b) show that after 3 Gyr the cold 
gas settles down in form of a stable torus, with no further 
condensation of extended cold gas. This is inconsistent 
with observations which show that about a third of cool- 
core clusters show Ha filaments extending out to 10s of 
kpc from the center (McDonald et al. 2010). The bot¬ 
tom panels in Figure 3 from our fiducial run show that 
the torus is unsteady even at late times with extended 
cold gas condensing out till the end of our run. We com¬ 


pare our results in detail with Li & Bryan (2014b) in 
section 4.1. 

3.1.4. Velocity and space distribution of cold gas 

We find it very instructive to classify the cold gas into 
two components: most of the mass is in the rotationally- 
dominant gas at < 5 kpc; a smaller fraction is in 
a radially dominant component spread over 20 kpc. 
Figure 4 shows the velocity and space distribution of 
rotationally (the left two panels; d 2 M/d\n \V(f,\dr and 
d 2 M/dln\v r \dr-. \v$\ > |w r |) and radially (the right two 
panels; d 2 M/d\n\v t f,\dr and d 2 M/d\n\v r \dr\ lu^l < |u r |) 
dominant cold (T < 5 x 10 4 K) gas, averaged from 1 
to 4 Gyr. The rotationally dominant gas distribution 
(lu^l > |v r |; two left panels in Fig. 4) shows two peaks 
at ±150 — 200 km s _1 and r « 1.5 — 4 kpc, cor¬ 
responding to the cold tori seen in Figure 3. The radial 
velocity is < 100 km s _1 . 

The distribution of the radially dominant cold gas in 
Figure 4 is quite different from the rotationally dominant 
gas. In addition to a larger radial extent, the radial ve¬ 
locity of the radially dominant component is much larger, 
going up to ±600 km s — , much larger than the maxi¬ 
mum azimuthal speed. The radial velocity of the closer 
in gas (< 5 kpc) is even larger for the outflowing (iy > 0) 
component because it is dredged up by the fast jet mate¬ 
rial; tiny mass in the cold gas is seen to reach a velocity 
close to Uj e t = 0.1c. The mass in the in-falling radially- 
dominant cold gas is ~ twice that of the outgoing cold 
gas. 

Figure 5 shows the 1-D velocity distribution of the cold 
gas averaged from 1 to 4 Gyr. The two large, sharp 
peaks correspond to the massive clockwise rotating cold 
torus. The radially-dominant component (|iv| > |u^|) 
shows a prominent high velocity tail in the positive di¬ 
rection. The negative velocity component for velocities 
larger than 300 km s _1 is also dominated by the radially 
in-falling (rather than rotationally dominant) gas, some¬ 
times affected by the fast jet back-flows. The maximum 
velocity peak of the radially and rotationally dominant 
cold gas coincide at ss 150 — 200 km s _1 , corresponding 
to the circular velocity at ~ 5 kpc. 

Figure 6 shows the relationship of in-falling and outgo¬ 
ing cold gas at small scales (5 kpc) and AGN jet activity. 
The cold outflow rate shows large spikes coincident with 
a sudden rise in AGN jet power, implying that cold gas 
observed with large velocity (inf Figs. 4 & 5) is dredged 
up by fast moving jets. The coincidence is particularly 
strong when a massive cold gas torus is present at small 
scales. For steady cooling in absence of angular momen¬ 
tum, we expect the mass inflow rate at 5 kpc and 1 kpc 
to closely follow each other. This is, however, not the 
case (especially around 2 Gyr) when majority of in-falling 
cold gas crossing 5 kpc is incorporated in the rotating 
cold torus, instead of accreting through 1 kpc. Also note 
that the outflowing cold gas is in form of very short-lived 
massive spikes, but the inflowing cold gas is smoother. 
The interpretation is that the outflowing cold gas is as¬ 
sociated with the AGN-uplifted cold torus gas, and the 
infalling cold gas is because of local thermal instability 
in a gravitational field. 

3.1.5. Cooling & heating cycles 
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Fig. 3.— The 2-D (z = 0) contour plots of electron number density (in cm -3 ) in the mid-plane of the very inner region at different times 
for the fiducial 3-D run, with the projection of the velocity unit-vector represented by arrows. The top-left panel shows the beginning of 
the infall of cold gas with random angular momentum. The top-second left panel shows an anti-clockwise transient torus. All times after 
this show a clockwise torus in the mid-plane which waxes and wanes because of cooling and AGN heating cycles. Even at late times the 
cold torus is not stable and gets disrupted by jets. 
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Fig. 4.- The velocity-radius distribution of the cold gas (T < 5 x 10 4 K) mass averaged from 1 to 4 Gyr. The top-left panel shows the 
v^ — r mass distribution ( dln \ v ^\ d in~ > ^ v 4> = A v r = 20 km s —1 , Ar = 0.5 kpc are the bin-sizes) for the rotationally-dominant (|u^>| > |u r |) 
gas and the bottom-left panel shows the v r — r distribution for the same gas. The top-right panel shows the — r distribution for the 
radially-dominant (|iy| > \v^\) cold gas and the bottom-right panel shows the v r — r ( ^i n ^i nr ) distribution for the same gas. Some 
of the salient features are: the rotationally-dominant cold gas, which is concentrated mainly within 5 kpc, is more abundant by a factor 
~ 40 than the radially-dominant gas; the dominant rotationally-supported clockwise cold torus with a negligible radial velocity (see Fig. 
3) is clearly visible in the two left panels; the radially-dominant cold gas (with \v r \ > 1^1) is much more radially extended, going out to 
25 kpc; the bottom-right panel shows that the in-falling ( v r < 0) cold gas dominates over the outgoing cold gas (by a factor « 2) and that 
the outgoing cold gas at < 5 kpc extends to very large velocities. 
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Fig. 5.— The velocity distribution of cold gas for the 3-D fidu¬ 
cial run with respect to the radial and azimuthal velocities. Also 
shown are the rotationally (|v^| > |i; r |) and radially (|i; r | > |tty|) 
dominant components. At large velocities the total and radi¬ 
ally/rotationally dominant components coincide but at small ve¬ 
locities they do not, as expected (at low velocities the other compo¬ 
nent of velocity dominates the mass budget). Also shown (dashed 
line) is the radial velocity distribution for the 2-D fiducial run; the 
azimuthal velocity is zero for 2-D axisymmetric runs. 


10000 


.1000 - 


g 100 ; 


T3 

0) 

N 


10 


0.1 


•^in,cold,5kpc(-^©y r 
•^out,cold,5kpc(-^©y r 1) 


Macc(M©yr - 1 ) 1 

i-tepowdr (111 41 nrg s ') 



1.5 2 

time(Gyr) 


Fig. 6.— The mass inflow (green short-dashed line) and outflow 
(red solid line) rates in the cold phase measured at 5 kpc as a 
function of time. Also shown are the jet power (normalized to 
10 41 erg s _1 ; how jet power is calculated is described in section 
3.1.5) and the mass accretion rate at 1 kpc (Eq. 6). Note that 
the largest spikes in the cold outflow rates are mostly associated 
with a sudden rise in jet energy, indicating cold gas uplifted by 
jets. The difference between the mass inflow rate at 5 kpc and at 
1 kpc (clearly noticeable at ~ 2 Gyr) leads to the build-up of the 
rotating cold torus (see Fig. 3). Also note that extended cold gas 
and jets are absent at ~ 3.75 Gyr. 


One of the distinct features of the cold feedback 
paradigm is that we expect correlations in jet power, 
cold gas mass, mass accretion rate, min(t coo i/fff), core 
entropy, etc. The observations indeed show such corre¬ 
lations (e.g., Figs. 1, 2 in Cavagnolo et al. 2008; see 
also Voit & Donahue 2014; Sun 2009; McDonald et 
al. 2011a). In Figure 7 we make ‘phase-space’ plots 
of jet power and cold-gas mass (total and the radially- 
dominant component) as a function of min(f coo i/tff) for 
our fiducial 3-D run. 


Evaluating min(t coo \/tg): The ratio t coo \/tg is calcu¬ 
lated by making radial profiles of emissivity-weighted 
(only including plasma in the range of 0.5 to 8 keV) in¬ 
ternal energy and mass densities. They are combined to 
calculate f C ooi = (3/2)n/csT/(n e niA), and t coo \/tg is cal¬ 
culated by taking its ratio with the free-fall time based 
on the NFW potential (Eq. 4; iff = pr/g] 1 / 2 where 
g = d$>/dr). The broad local minimum in t coo \/tg pro¬ 
file is searched going in from the outer radius and is used 
as min(t cool /t ff ). 

Evaluating jet power: The jet power is also calculated 
in a novel way, which is close to what is done in obser¬ 
vations. 3 We consider the grids with mass density lower 
than a threshold value (chosen to be 0.17 times the initial 
minimum density in the computational volume; results 
are insensitive to the exact value of the threshold den¬ 
sity) to belong to the jet/bubble material, and we simply 
volume-integrate the internal energy density of all such 
cells to calculate the jet energy (only considering thermal 
energy; we use 7 = 5/3 for the jet material because it is a 
non-relativistic hot gas in our simulations; in reality, rel¬ 
ativistic particles with 7 = 4/3 are a major component 
of jets). This density-based definition of the jet material 
coincides with the visual appearance of the jet. The jet 
energy is divided by an estimate of the bubble lifetime 
(chosen to be 30 Myr, of order the dynamical/buoyancy 
time at 10 kpc; see Table 3 in Birzan et al. 2004) to 
arrive at the jet power. For simplicity, we use the same 
value of the bubble lifetime at all times in all our runs. A 
trivial definition of jet power, in which it is proportional 
to the instantaneous accretion rate at 1 kpc, is given by 
Eq. 6 as 3.4 x 10 42 M acc ( Mgyr -1 ) erg s _1 . Assuming 
this conversion, Figure 6 shows that the two estimates of 
jet power are comparable in magnitude but vary rather 
differently with time. This is because, while M acc is an 
instantaneous quantity varying on a dynamical timescale, 
our jet power is based on the jet thermal energy which 
is an integrated quantity. 

We anticipate cycles in the evolution of min(t CO oi/tff) 
and jet power or the radially-dominant cold gas. Imag¬ 
ine that there is no accreting cold gas at the center; in 
this state without heating the core is expected to cool 
below t C ooi/iff~ 10 (because accretion rate in the hot 
mode is small). The t coo \/tg< 10 state is prone to cold- 
gas condensation and enhanced feedback heating if e is 
sufficiently high. Energy injection leads to overheating 
of the core and an increase in t coo \/tg] since condensa¬ 
tion/accretion is suppressed in this state, both jet power 
(because of adiabatic/drag losses) and radially-dominant 
cold gas mass are reduced in this state of t coo i/tg> 10. 
Eventually the core cools again and the cycle starts 
afresh. 

The left panel of Figure 7 shows one of the many 
jet cycles in our fiducial 3-D cluster run. On average 
jet power vs. min(t coo i/tff) evolves in form of clock¬ 
wise cycles of various widths (a measure of the range of 
min[f coo i/tfj] before and after the jet event) and heights 
(jet power). Generally, a smaller t coo \/tg leads to a larger 

3 Observers calculate the bubble/cavity mechanical power by 
assuming it to be in pressure balance with the background ICM 
and by using a size and an age estimate for the bubble (e.g., see 
Birzan et al. 2004). Indeed, our bubbles are in pressure balance 
with the ICM, as seen in Fig. 1. 
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Fig. 7.— The variation of jet power, total cold gas mass and the radially-dominant cold gas mass as a function of min(t coo l/^ff) for the 
fiducial 3-D run from 2.2 to 3.2 Gyr. Color shows the evolution of cluster in time. While jet power (left panel) and radially-dominant cold 
gas mass (right panel) show clockwise cycles with min(t coo i/tff), the total cold gas mass (middle panel) simply builds up in time. Notice 
the linear scale for the total cold gas mass instead of a log scale for the other two cases. 


mass accretion rate and a larger jet energy, and therefore 
larger overheating and a larger min(t coo i/tff). Since the 
efficiency of our fiducial run is rather small (e = 6 x 10 -5 ), 
the cluster core remains with t coo i/tff< 20 at most times. 
In section 3.2.2 we discuss the dependence of our results 
on jet efficiency (e). 

The middle panel of Figure 7 shows the total mass in 
cold gas (most of which is in the cold rotating torus) 
as a function of min(t coo i/tfj). We see the mass in the 
cold torus building up in time. We can easily see that 
the total cold gas mass simply builds up in time (see the 
green dashed line in the upper panel of Fig. 9), and 
is uncorrelated with min(f coo i/tff). The right panel of 
Figure 7 shows the mass in the radially-dominant cold 
gas (with \v r \ > \vj,\) as a function of min(f coo i/fff). This 
panel also shows clockwise cycle like jet power shown in 
the left panel. A larger radially-dominant cold gas mass 
generally implies a higher accretion rate and a larger jet 
power, but the features in jet and cold gas cycle are not 
always varying in an identical fashion. While the global 
evolution in phase space is clockwise, there is haphazard 
evolution at smaller timescales (e.g., between 2.5 to 2.9 
Gyr). 

3.2. The 2-D runs 

The 3-D simulations are very expensive compared to 
the 2-D ones, not only because the number of grid cells 
is larger but also because the CFL time step is much 
smaller. The CFL time step in 3-D is dominated by cells 
close to the polar regions (9 = 0, 7r) and oc r sin 9A<p ss 
rA9A<p/2, much smaller than in 2-D (oc rA9). Our 
256 x 128 x 32 (3-D) runs have 8 times more grid cells 
compared to our 512 x 256 (2-D) runs and the CFL time 
step is ss 0.2 times smaller, making our 3-D runs 40 times 
more expensive than the 2-D ones. Therefore, for scans 
in various parameters (halo mass, jet efficiency, etc.), 


only 2-D axisymmetric simulations are practical. How¬ 
ever, the key drawback is that the initially non-rotating 
gas cannot gain angular momentum in axisymmetry, and 
thus 2-D simulations do not show the formation of a rota- 
tionally supported torus. But, as we discuss shortly, the 
suppression of cooling flow, nature of radially-dominant 
cold gas, etc. are very similar in 2-D and 3-D. 

In this section we describe different variations on the 
fiducial setup for our 2-D simulations. Section 3.2.1 com¬ 
pares the results from 2-D and 3-D simulations with cool¬ 
ing and AGN feedback. Section 3.2.2 studies the effect 
jet feedback efficiency and the halo mass on the proper¬ 
ties of jets and cold gas. 

3.2.1. Comparison with 3-D 

Since 3-D simulations are substantially more time con¬ 
suming compared to the 2-D axisymmetric ones, it will 
be very useful if some robust inferences can be drawn 
from these faster 2-D computations. To compare the 3-D 
simulations with their 2-D counterparts, we have carried 
out the fiducial 3-D simulation in 2-D with identical pa¬ 
rameters (the initial density perturbations in 2-D runs 
are the same as the perturbation for the <j> = 0 plane in 
3-D). 

Figure 5 compares the time-averaged velocity distribu¬ 
tion of cold gas in 2-D and 3-D simulations. Since the 
azimuthal velocity vanishes in 2-D axisymmetric simu¬ 
lations, we compare the radial velocity distribution of 
cold gas in 2-D simulations with the radially-dominant 
(|w r | > \v$\) component in 3-D. While the outflowing 
cold gas has a similar distribution in 2-D and 3-D, the 
inflowing gas is more dominant in 2-D relative to 3-D be¬ 
cause in 3-D a lot of this in-falling cold gas slows down 
and becomes a part of the rotating cold torus. 

Table 1 shows that the mass accretion rate through 
the inner boundary for the fiducial runs in 2-D and 3- 
D are comparable. Unlike in 3-D, we note that there is 













12 


D. Prasad, P. Sharrna, A. Babul 



Fig. 8. — The mass accretion rate relative to the cooling flow 
value as a function of jet efficiency (Eq. 6) for cluster and massive 
cluster runs (both 2-D and 3-D). The accretion rate is suppressed 
more for a lower halo mass at a fixed e. 


substantial cold gas sticking to the poles in 2-D due to 
numerical reasons. Similarly, in 3-D there is a physical 
accumulation of cold gas in form of a rotating torus. 

Figure 8 shows the average mass accretion rate through 
the inner radius of our simulation volume (M acc ) relative 
to the cooling flow rate (M c f). The suppression factor 
(Macc/M c {) for 3-D cluster simulations (with e = 6 x 
10~ 5 , 5 x 1CT 4 , 0.01) is comparable to 2-D. 

Figure 9 shows various important quantities, such as 
jet power, cold gas mass, mass accretion rate through the 
inner boundary, as a function of time for the fiducial 3-D 
(upper panel) and 2-D (bottom panel) runs. Encourag¬ 
ingly, various quantities, except the total cold gas mass, 
show similar trends with time in 2-D and 3-D. The to¬ 
tal cold gas mass is much larger in 3-D because of the 
formation of a massive cold torus which is absent in ax- 
isymmetry. 

In both 2-D and 3-D runs min(f coo i/iff) varies in the 
range 1 to 10, and is roughly anti-correlated with M acc 
and jet power. The maximum jet power goes up to 
~ 10 45 erg s _1 in both cases. The mass accretion rate 
and hence feedback power injection (Eq. 6) is more spiky 
in 2-D (can go above 100 M 0 yr _1 for some times) be¬ 
cause, unlike in 3-D, the cold gas that is accreted in 2-D 
covers full angle 2tt in <f> because of axisymmetry. The 
jet power, which is calculated by measuring the instan¬ 
taneous jet thermal energy, depends on the average mass 
accretion rate over < 100 Myr rather than the instanta¬ 
neous value. Another difference between 2-D and 3-D is 
that cold gas can be totally removed (through the inner 
boundary) after strong feedback jet events in 2-D but 
this never happens in 3-D; cold gas (even the radially- 
dominant component) is present at all times because it 
is very difficult to evaporate/accrete the massive rotat¬ 
ing cold torus. There definitely is a depletion in the 
amount of radially-dominant cold gas after a strong feed¬ 
back event in 3-D (at ~ 3.5 Gyr in the top panel of Fig. 

9 )- 

Although we have not explicitly shown jet power and 
cold gas ‘phase-space’ plots for our 2-D fiducial run, we 
have verified that it shows cycles similar to the 3-D run 
(left & right panels of Fig. 7). Indeed, Figure 9 in¬ 
dicates that the 2-D runs should also show clock-wise 
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Fig. 9. — Various quantities (jet power, cold gas mass, radially- 
dominant cold gas mass, min[t coo i/tff], M acc ) as a function of time 
in the fiducial 3-D (top panel) and 2-D (bottom panel) cluster 
simulations. The data are sampled every 10 Myr. The cycle shown 
in Figure 7 are based on the top panel. All quantities, except total 
cold gas mass, are statistically similar in 3-D and 2-D. The total 
cold gas mass, which is dominated by the cold torus, is much larger 
in 3-D and builds up in time. However, the radially-dominant cold 
gas (| Ur | > |u^|) mass in 3-D is similar to the total cold gas mass 
in 2-D. 


cycles in jet energy and cold gas mass as a function of 
min(f coo i/tff). These cycles just reflect the sudden rise 
in the accretion rate (M acc ) and jet power due to cold 
gas condensation and slow relaxation to equilibrium af¬ 
ter overheating (notice the fast rise and slow decline in 
jet energy for individual jet events in both panels of Fig. 
9). 


3.2.2. Dependence on jet efficiency & halo mass 

Till now we have discussed the fiducial cluster simula¬ 
tion with a small feedback efficiency e = 6 x 10 -5 . In 
this section we study the influence of jet efficiency (e) 
and halo mass (M 2 oo) on various properties of the clus¬ 
ter core. Overall, we find that the effect of an increasing 
halo mass is similar to that of a decreasing feedback ef¬ 
ficiency. We compare the efficiencies (e) ranging from 
6 x 10~ 6 to 0.01. We consider two halo masses: a cluster 
with M 200 = 7 x 10 14 Mq and a massive cluster with 
M 200 = 1.8 x 10 15 M 0 . 

Table 1 and Figure 8 show that the feedback efficiency 
ofe = 6xl0 _5 is able to suppress the cooling flow by 
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cluster, e = 6 x 10 -6 ■•■»■■■ massive cluster, e = 6 x 10 -5 

cluster, e = 6 x 10 -5 '■ massive cluster, e = 5 X 10 -4 
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) 05 1 t5 2 2T5 3 t5 \ A5 5 


time(Gyr) 

Fig. 10.- The mass accretion rate through the inner radius 
(smoothed over 50 Myr) as a function of time for different 2-D 
runs. A lower efficiency and a more massive halo lead to a larger 
accretion rate. 

about a factor of 10 for a cluster but only by a factor 
of 4 for a massive cluster. This implies that a larger ef¬ 
ficiency is required to suppress a cooling flow in a more 
massive halo. We note that the pure cooling flow accre¬ 
tion rate decreases with a decreasing halo mass because 
of a smaller amount of gas in lower mass halos (see the 
values enclosed in brackets in Table 1). 

Figure 8 shows that the suppression factor (M a cc/M c f) 
is smaller for the massive cluster, and scales as e -2 / 3 for 
both cluster and massive cluster runs (see also Table 1). 
A decrease in the accretion rate with an increasing e is 
not a surprise; a higher feedback efficiency heats the core 
more and maintains t coo \/tg> 10 at most times, resulting 
in only a few cooling/feedback events. While the average 
jet power (~ eM acc c 2 oc e 1 / 3 ) increases with an increasing 
e, the core X-ray luminosity decreases. This implies that 
feedback heating and cooling do not balance each other 
at all times. Heating dominates cooling just after jet 
outbursts and cooling dominates in absence of infalling 
cold gas when t coo \/ts slowly decreases from a value 
> 10. Thus, for a larger e, for which a cluster spends 
more time in a hot/dilute state, the X-ray emission from 
the core is expected to be smaller (c.f., Fig. 12). 

Figure 10 shows the mass accretion rate (averaged over 
50 Myr bins) as a function of time for our 2-D cluster and 
massive cluster runs with different feedback efficiencies. 
The solid red line corresponds to the fiducial 2-D cluster 
run (with e = 6x 10 -5 ). The green dotted line with a 
marker, which corresponds to a ten times lower efficiency 
(e = 6 x 10 -6 ), shows an accretion rate comparable to a 
cooling flow at most times (see also Table 1). The clus¬ 
ter run with ten times higher efficiency (e = 5 x 10 -4 ), 
indicated by black dot-dashed line, shows an average ac¬ 
cretion rate of 5.1 M 0 yr _1 (about a fifth of the fiducial 
2-D run; see Table 1); there are far fewer spikes in M acc 
compared to the fiducial run. Similar trends are observed 
for the massive cluster runs with e = 6 x 10 -5 (magenta 
dotted line) and 5 x 10~ 4 (blue double-dot-dashed line). 
The number of M spikes in Figure 10 are smaller for 
lower halo mass and higher feedback efficiency because 
of larger overheating and a longer recovery time after a 


precipitation-induced jet event. 

Figure 11 shows the time averaged (from 4 to 5 
Gyr) and emissivity weighted (0.5-8 keV) 1-D profiles 
of several key quantities for 2-D cluster and massive 
cluster runs with different efficiencies: entropy (K = 

Tkev/fie^ 3 ), tcooi/tff , number density and temperature. 
All profiles look similar to what is seen in observations. 
The entropy profile flattens toward the center but the en¬ 
tropy core is prominent only for the higher efficiency (e) 
runs; a ‘core’ with a constant t coo \/ts is more prominent 
for lower e and for the massive cluster. As expected, the 
density is lower and the temperature is higher for a larger 
feedback heating efficiency. For all efficiencies tempera¬ 
ture increases with the radius (except for e = 5 x 10~ 4 
which is almost isothermal just after a jet outburst; see 
the top-right panel of Fig. 12), as seen in the observa¬ 
tions of cool-core clusters. 

Compared to the cluster runs, the entropy for the mas¬ 
sive cluster is higher at larger radii in Figure 11 be¬ 
cause the initial entropy was scaled with the halo mass 

(K oc M200 ; see section 2.2). The entropy profiles for the 
massive cluster runs for the two efficiencies are similar; 
entropy keeps on decreasing as we go toward the center 
(forming a ‘core’ in t coo i/tg), more so for e = 6 x 10~ 5 . 
As we saw with the mass accretion rate in Figure 10, 
the effect of increasing the efficiency is similar to that of 
decreasing the halo mass. This is expected, as the mass 
accretion rate for lower mass halos is smaller, and the 
increase in jet efficiency and the consequent higher jet 
power suppresses accretion. 4 

Another point to note in Figure 11 is that the pro¬ 
files are rather similar for the massive cluster runs with 
e = 6 x 10~ 5 and e = 5x 10 -4 . The bottom panels of 
Figure 12 show that jet events between 4 to 5 Gyr are not 
able to raise min(t coo i/iff) much above 10 for these cases. 
However, top panels of Figure 12 and the bottom panel 
of Figure 9 show that between 4 to 5 Gyr min(f coo i/fff) 
increases with an increasing e. Therefore, the core en¬ 
tropy (density) for the cluster runs increases (decreases) 
with an increasing e. Note that the core entropy for the 
cluster runs with a larger efficiency are not always higher; 
its only when the core is in the part of the heating cycle 
with min(t coo i/tff)> 10. 

Figure 12 shows various quantities (jet power, cold gas 
mass and min[£ coo i/ffj]) as a function of time for 2-D clus¬ 
ter (e = 10“ 4 and 5 x 10 -4 ; also see the 2-D cluster run 
with e = 6 x 10 -5 in the bottom panel of Fig. 9) and 
massive cluster (e = 6 x 10 -5 and 5 x 10 -4 ) runs. The 
first point to note is that the number of jet events (and 
hence the number of cycles; e.g., see Fig. 7) is smaller for 
a higher efficiency and a lower halo mass. Another is that 
the peaks in jet power and min(t coo i/tfj) for a higher effi¬ 
ciency are larger, resulting in overheating and longer du¬ 
rations for which cold gas and jet power are suppressed. 
Stronger overheating after jets in higher efficiency (and 
lower halo mass) runs results because, while the number 
of cold accretion events are smaller (compared to lower 
efficiency or a larger halo mass), the mass accretion rate 
during the multiphase cooling phase is similar (see Fig. 
10), generally giving larger heating (Eq. 6). 

As with the mass accretion rate (see the spikes in Fig. 

4 We thank the referee for the suggestion to highlight this point. 
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Fig. 11. — Emissivity-weighted (considering plasma in the range 0.5 — 8 keV), time (averaged between 4 and 5 Gyr) and angle-averaged 

profiles for 2-D runs as a function of radius scaled to the outer radius (r max ): entropy (K = T^ e y/n^ 3 top-left panel); t coo l/ffF (top-right 
panel); electron density (n e ; bottom-left panel); and temperature (in KeV; bottom-right panel). Both min(t coo [/tfj) and core entropy 
decrease for a lower efficiency or a larger halo mass. Temperature is higher for a higher efficiency, but a cool core (a temperature increasing 
with radius) is preserved for all cases, except the highest efficiency run. Spikes in the cluster run with e = 5 X 10“ 4 signify that there is 
cool, low entropy gas present in the core from 4 to 5 Gyr. 

10), for a fixed e the number of cooling/jet events are 
larger for a massive cluster. While t coo i/ts is < 10 and 
cold gas is present at most times for the massive clus¬ 
ter run with e = 6 x 10” 5 , there are longer periods with 
min[t coo i/fff]> 10 and lack of cold gas for the cluster run 
(see bottom panel of Fig. 9). The jet events are more 
disruptive (as measured by the rise in min[t coo i/fff] af¬ 
ter a jet event) in the lower mass halo because the jet 
power is relatively large but the hot gas mass is smaller 
(compare the right panels of Fig. 12). 

4. DISCUSSION & ASTROPHYSICAL IMPLICATIONS 

The cold mode accretion model,' 1 in which local ther¬ 
mal instability leads to the condensation and precipita¬ 
tion of cold gas and enhanced accretion on to the SMBH, 
has emerged as a useful framework to interpret various 
properties in cores of elliptical galaxies, groups, and clus¬ 
ters (e.g., Pizzolato & Soker 2005; Sharrna et al. 2012a; 

Gaspari et al. 2012; Li & Bryan 2014b). However, there 
are several unresolved problems: e.g., the role of angular 

5 Here we use the label “cold mode accretion” to refer to the 
capture and accretion of cold clouds by SMBH, and not the cos¬ 
mological accretion of cold gas sometimes invoked in halos less 
massive than 10 12 Mq (Birnboim & Dekel 2003). 


momentum transport, self-gravity and cloud-cloud col¬ 
lisions in accretion on to the SMBH (e.g., Pizzolato & 
Soker 2010; Hobbs et al. 2011; Babul et al. 2013; 
Gaspari et al. 2014); relative contribution of cold gas 
at ~ 1 kpc to SMBH accretion and star-formation; the 
role of thermal conduction in thermal balance and cold- 
gas precipitation (Wagh et al. 2014; Voit & Donahue 
2014); the exact mechanism (turbulent mixing, weak 
shocks; e.g., see Fabian et al. 2003; Dennis & Chandran 
2005; Banerjee & Sharrna 2014) via which the mechan¬ 
ical power of the jet/cavity is dissipated and distributed 
throughout the core; and the scaling of various processes 
with the halo mass. 

The key feature of cold mode accretion, unlike the hot 
mode, is that the mass accretion rate increases abruptly 
as t coo \/tff becomes smaller than a critical value close 
to 10. This leads to a strong feedback heating, which 
temporarily overheats the cluster core. The hot mode 
feedback, in form of Bondi accretion onto the SMBH, on 
the other hand, is not an abrupt switch and increases 
smoothly with an increasing (decreasing) core density 
(temperature). In section 4.1 we discuss the success 
of the cold accretion model and compare with previous 
simulations. In section 4.2 we compare with the recent 











Cold gas and AGN jet feedback in cluster cores 


15 



time (Gyr) time (Gyr) 


Fig. 12.— Jet energy, cold gas mass, and min(t coo i/tff) as a function of time for 2-D runs with different efficiencies (e = 10 4 , 5 X 
10 4 , 6 X 10 5 ) and halo masses (M 200 = 7 X 10 14 , 1.8 X 10 15 Mq). Note that the jet energy and cold gas mass are scaled differently in 
different panels. A smaller efficiency or a larger halo mass leads to many accretion and jet feedback cycles. 


exquisite cold-gas observations and with statistical anal¬ 
yses of X-ray and radio observations of cluster cores. 

4.1. Comparison with previous simulations 

There are two broad categories of jet implementations 
described in the literature: first, where the jet mass, mo¬ 
mentum and energy are injected via source terms (e.g., 
Omma et al. 2004; Cattaneo & Teyssier 2007; Li & 
Bryan 2014a; Gaspari et al. 2012); second, where mass 
and energy are injected as flux through an inner bound¬ 
ary (e.g., Vernaleo & Reynolds 2006; Sternberg, Pizzo¬ 
lato, & Soker 2007). We use the former approach, which 
has generally been more successful. In this approach, the 
sudden injection of kinetic energy after cold gas precipi¬ 
tation leads to a shock which not only expands vertically 
but also laterally, perpendicular to the direction of mo¬ 
mentum injection. This lateral spread of jet energy and 
vorticity generation due to interaction with cold clumps 
help in coupling the jet energy with the equatorial ICM. 
In the flux-driven approach the jet pressure is usually 
taken to be the same as ICM pressure and the jet drills 
a cavity without expanding laterally in the core. Thus, 
coupling of the jet power is not very effective, unless the 
jet angle is very broad (Sternberg, Pizzolato, & Soker 
2007). 

Our jet modeling is similar to the earliest works such 
as Omma et al. (2004); Omma & Binney (2004), which 
inject jet mass, momentum and kinetic energy via source 
terms. However, this work focussed on the effect of a sin¬ 
gle jet outburst with a fixed power and did not include 
cooling; the simulations were run for short times (< 500 
Myr). Cattaneo & Teyssier (2007) also implement jets 


using kinetic and thermal energy injection, and run for 
cosmological timescales. However, they use Bondi pre¬ 
scription for accretion on to the SMBH, and hence their 
jet power input is tied to cooling and accretion at the 
center. Since the Bondi radius cannot be resolved in 
their simulations, they compute the Bondi accretion rate 
based on the density and temperature at a larger radius. 
Sometimes (e.g., in Dubois et al. 2010) the Bondi accre¬ 
tion rate evaluated at large radii is artificially enhanced 
by a factor ~ 100 in order to match feedback heating with 
cooling. Bondi accretion is only applicable for a smooth, 
non-rotating gas distribution, and not for clumpy multi¬ 
phase gas which can accrete at a much higher rate (e.g., 
Gaspari et al. 2013; Sharma et al. 2012a). 

Cielo et al. (2014) have studied the detailed structure 
and thermodynamics of source-term driven cylindrical 
jets, of different densities and temperatures, interacting 
with the ICM but they run for less than 10 Myr. Like 
us, they also highlight the importance of hot back flows 
in regulating the central ICM. 

Another set of simulations inflate cavities using jets 
driven by fluxes of mass and momentum at the inner ra¬ 
dial boundary (rather than using source terms like us; 
in cluster context, see Vernaleo & Reynolds 2006; Stern¬ 
berg, Pizzolato, & Soker 2007; for MHD modeling of 
the Crab nebula jet, see Mignone et al. 2013). Vernaleo 
& Reynolds (2006) injected momentum (and kinetic en¬ 
ergy) via the inner radial boundary, with an opening an¬ 
gle of 15°, in form of 100 times hotter gas but in pressure 
equilibrium with the ICM. Their jets just drill through a 
narrow channel without coupling to the catastrophically 
cooling core. 
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Sternberg, Pizzolato, & Soker (2007) advocated wide 
(with opening angle > 50°) boundary-driven jets, such 
that the jet is not as fast, and can lead to vortices and 
substantial mixing in cluster cores. However, since their 
simulations are not run for many cooling times, its un¬ 
clear if wide jets and can indeed balance cooling for cos¬ 
mological times. Moreover, the fat jets may not repro¬ 
duce the observed morphologies of thin jets and fat bub¬ 
bles. Using the boundary injection approach, Heinz et 
al. (2006) emphasize the importance of the dynamic 
ICM in redistributing jet energy but they also run for 
less than a cooling time. 

Recent numerical simulations of AGN-driven jets (Gas- 
pari et al. 2012; Li & Bryan 2014a,b) have been quite 
successful in producing several observed features such as, 
the lack of plasma cooling below a third of the ICM tem¬ 
perature (Fig. 11 in Li & Bryan 2014b), suppression 
of cooling and accretion in the core (by a factor of 10- 
100 relative to a cooling flow), maintenance of cool-core 
structure even with strong intermittent jet events, forma¬ 
tion of an angular momentum supported cold-gas torus, 
viability of AGN feedback from elliptical galaxies to mas¬ 
sive clusters. Our simulations are different from these 
recent works, which use mesh refinement in a cartesian 
geometry, in that we use a spherical coordinate system. 
We have also tried to push the AGN feedback efficiency 
toward the lower limit which is still able to suppress a 
cooling flow. We find that an efficiency e = 6xl0 _5 is 
able to suppress a cluster cooling flow by a factor of 10. 

Like us, Gaspari et al. (2012) and Li & Bryan (2014b) 
also make an estimate of the mass accretion rate on to 
the SMBH. Gaspari et al. (2012) consider a spherical 
shell of radius 0.5 kpc and calculate the mass accretion 
rate (M acc ) due to infalling gas. Li & Bryan (2014b); Li 
et al. (2015) calculate the mass accretion rate (M ac c) by 
dividing the cold gas mass within 0.5 kpc by 5 Myr (of 
order the dynamical timescale). Our estimate of M acc is 
similar to Gaspari et al. (2012), except that we calculate 
it at 1 kpc. Only a small fraction of M acc is expected to 
be accreted onto the SMBH; thus, the efficiency factor (e) 
in Eq. 6 takes into account both the fraction of M acc that 
is accreted by the SMBH and the efficiency of converting 
SMBH accretion into jet mechanical energy. 

While our jet feedback implementation is very similar 
to Gaspari et al. (2012), our results differ in some key re¬ 
spects. The main difference is that we see extended cold 
gas and jet/cold-gas cycles even at late times (see Figs. 
1, 3, 7, 9). Like Li & Bryan (2014b), in Gaspari et al. 
(2012) there is a long-lived rotationally supported torus 
at few kpc and the extended multiphase gas is lacking at 
later times (see their Figs. 10 & 11). The main reason 
for the absence of extended cold gas and strong jets at 
late times in previous simulations is a large feedback ef¬ 
ficiency. A larger feedback efficiency leads to very strong 
feedback heating at early times, and the core reaches 
rough thermal balance in a state of t C ooiAff> 10 with no 
fresh extended (radially dominant) cold gas condensing 
at late times. Since a large fraction of cool core clus¬ 
ters show extended cold gas (McDonald et al. 2011a), a 
smaller value of feedback efficiency seems more consistent 
with observations. 

To solve the problem of a steady cold torus present at 


late times, very recently, Li et al. (2015) have incor¬ 
porated the depletion of cold gas via star formation in 
the core, but they adopt the same feedback prescription 
as in their previous works. Star formation exhausts the 
amount of cold gas within 0.5 kpc, suppresses AGN heat¬ 
ing, and leads to a cooling event after which cold gas con¬ 
denses again. Thus, they obtain three cooling-feedback 
cycles in their fiducial run. In their AGN feedback pre¬ 
scription, star formation efficiency primarily determines 
the frequency of cooling/heating cycles. In contrast, our 
cycles are determined by the AGN feedback efficiency 
and the halo mass (more cycles for a massive halo and a 
smaller feedback efficiency). While there is some ongo¬ 
ing star formation in cool cluster cores, Li et al. (2015) 
form 3 x 10 11 M 0 in stars over 6.5 Gyr for their fidu¬ 
cial run corresponding to the Perseus cluster, a signifi¬ 
cant fraction of the mass of the BCG (brightest cluster 
galaxy; 8 x 10 11 M 0 ; Lim et al. 2008). This does not 
agree with semi-analytic models which suggest that 80% 
of the stars of BCGs are assembled before z = 3 (i.e., 
only 2 x 10 10 M 0 are expected to form over the past 12 
Gyr; De Lucia & Blaizot 2007). Moreover, the current 
star formation rate even in the most extreme cool core 
clusters is typically < 10 M 0 yr -1 (Hicks, Mushotzky, 
& Donahue 2010; McDonald et al. 2011b); the average 
star formation rate of Li et al. (2015) would correspond 
to an unacceptably high value, ~ 46 M 0 yr . 

We can directly compare our results with Gaspari et 
al. (2012) as our feedback prescription is similar. They 
tried jet efficiency factors of e = 6 x 10 , 0.01. With 

a much higher accretion efficiency (e = 0.01) compared 
to ours (~ 10~ 4 — 10 -5 ), Gaspari et al. (2012) get a 
larger suppression factor (M acc /M c f ~ 1 — 2%) compared 
to us, as expected. Gaspari et al. (2012) use a halo 
mass M 2 oo ~ 10 15 M 0 . Figure 8 shows how our results 
compare with Gaspari et al. (2012). The suppression 
factor of a massive cluster (M 2 oo = 1.7 x 10 15 M 0 ) for 
our fiducial e = 6 x 10 -5 is 20%, larger than their work. 
Suppression factor in our massive cluster (cluster) run 
for e = 0.01, as seen in Figure 8, is 3% (0.8%), in rough 
agreement with the results of Gaspari et al. (2012). 

4.2. Comparison with observations 

Now that we have done some comparisons with previ¬ 
ous simulations, in this section we compare our results 
with observations. We note that the observational com¬ 
parison may not be perfect because our simulations lack 
some physical processes such as magnetic fields and ther¬ 
mal conduction. These effects will be considered later. 
Moreover, observations suffer from projections effects, 
and our cluster parameters do not span as broad a range 
(of halo masses, entropy profiles at large radii, etc.) as 
encountered in observations. 

One of the most commonly studied ICM property is 
its entropy profile (e.g., Cavagnolo et al. 2009). Figure 
11 shows the time averaged (4-5 Gyr), X-ray emissivity 
weighted profiles of t coo \/tg and entropy as a function 
of radius for our various runs. We see that an entropy 
core (with a prominent local minimum in t C ooiAff) is a 
good approximation for systems in which t coo \/tg > 10 
and there is no extended cold gas. This state is common 
for simulations with larger efficiency and lower halo mass. 
However, the halos with t coo \/tg< 10, in which there is 
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Fig. 13.— Various important quantities measured at the same time (jet energy, cold gas mass, core entropy, and min[t coo i/tff]) plotted 
against each other from our 2-D cluster runs with different efficiencies. The data is sampled every 10 Myr. There is a strong correlation 
between the core entropy and min(£ coo i/%), especially at larger values of min(£ coo i/£ff). There is also a positive correlation between Ko-jet 
energy and min(t coo i/£ff )-jet power. Larger efficiency runs lead to a larger value of min(£ coo i/£ff) and Kq. Notice that cold gas is absent if 
min(t coo i/t ff )> 30. 



Fig. 14.— Important quantities measured at the same time (jet energy, cold gas mass, core entropy, and min[t coo i/%]) plotted against 
each other in our fiducial 3-D cluster run. As in 2-D runs (see Fig. 13), there is a strong correlation between Ko-mm{t coo \/tff ), Kq— jet 
power, and min(£ coo i/iff )-jet energy. The cold gas mass is high and becomes almost constant at later times as seen in the top panel of Fig. 
9. The color-coding corresponds to time. 
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lot of currently condensing and infalling cold gas, are 
consistent with a ‘core’ in t coo \/tg and a decreasing en¬ 
tropy toward the cluster center, albeit with a shallower 
slope. This is consistent with recent reanalysis of core 
entropy profiles (Panagoulia, Fabian, & Sanders 2014), 
which suggests that a double power-law entropy profile, 
with a shallower entropy in the center, better describes 
the ICM core. It will be useful to compare the behav¬ 
ior of central entropy as a function of min(t coo i/tff); we 
expect entropy cores for min(f coo i/tff)> 10 and slowly 
increasing entropy profiles for min(f coo i/tff)< 5. 

Figures 13 and 14 show the correlation between vari¬ 
ous important quantities for our 2-D cluster runs (with 
e = 6 x 10 -5 , 10 -4 , 5 x 10 -4 ) and the 3-D fiducial 
run, respectively. Data points sampled every 10 Myr 
are shown. The core entropy (Kq) is obtained by using a 
least squares fit to the emissivity-weighted 1-D entropy 
profile of gas in 0.5-8 keV range. In both these figures the 
strongest correlation is between K 0 and min(t coo i/tfj), as 
expected, because both these quantities depend on den¬ 
sity and temperature in a similar way ( K oc T/n 2 ' 3 and 
tcooi/^ffcx T 1/,2 /n; see Eq. 35 in McCourt et al. 2012); 
the relation is not one-to-one because Kq is determined 
by entropy near the center and min(t coo i/tff) by the be¬ 
havior at the core radius (beyond which density decreases 
sharply). 

The spread in A'o — min(t coo i/tff) correlation is larger 
for a lower Kq (or equivalently, min[t coo i/tff]; this is also 
seen in observational data shown in Fig. 4 in Voit & 
Donahue 2014) because a core with constant entropy 
is not a good description in that case and the entropy 
decreases inward (see top-left panel of Fig. 11). 

The correlation between various quantities in Figures 

13 and 14 are not particularly strong because of the hys¬ 
teresis behavior of various quantities (e.g., jet power, 
radially dominant cold gas mass) with respect to the 
core properties (Fig. 7). Figures 13 & 14 show that, 
in general, the jet power increases for a larger Kq (or 
min[f coo i/tff]), particularly for a larger core entropy. This 
is because a large jet power overheats the cluster core and 
raises its entropy. Other quantities do not show as strong 
correlations in these plots; cold gas mass increases with 
a lower entropy or a shorter cooling time, but jet energy 
and cold gas mass show a large spread relative to each 
other (since cold gas leads to increase in jet energy, which 
in turn suppresses cold gas mass). 

The 3-D run shown in Figure 14 prominently shows 
the sign of the massive torus at late times. Apart from 
this, there are no major differences in 2-D and 3-D. Also 
note that cold gas is missing in e = 5 x 10~ 4 2-D cluster 
run for t coo \/tg> 20 (Fig. 13; the same is expected for 
the radially dominant cold gas in 3-D). This is consistent 
with the observations of Cavagnolo et al. (2008), who 
find that Ha luminosity is suppressed for a core entropy 
> 30 keV cm 2 (corresponding to min[f C ooi/te] of about 
20; see the top-right panel of Fig. 13 & Fig. 4 in Voit 
& Donahue 2014). The onset of star formation in clus¬ 
ter cores also happens sharply below the same entropy 
threshold (Rafferty, McNamara, & Nulsen 2008). Figure 

14 shows that the core entropy and min(f coo i/tff) remain 
below 20 even if the instantaneous jet power is as high 
as 10 45 erg s _1 ; core entropy can be much higher (up to 
100 keV cm 2 ) for higher efficiency (see e = 5 x 10 -4 in 


Fig. 13). 

While Figures 13 & 14 show correlations between var¬ 
ious quantities at a given time, we are also interested 
in understanding causal relationships between various 
quantities such as cold gas mass, min(f coo i/fg), jet en¬ 
ergy. Figure 15 shows the temporal cross-covariance be¬ 
tween various important parameters (we take log before 
calculating cross-covariance). Various 2-D and 3-D sim¬ 
ulations with different efficiencies are plotted together. 
The trends which are common to all simulations are likely 
to be robust. Robust correlations among various quan¬ 
tities occur only with a time lag < 1 Gyr (typical core 
cooling/heating time). 

Top three panels in Figure 15 show that the cross¬ 
covariance between min(t coo i/tgf) decreases ~ 0.5 Gyr 
before zero lag and rises for 0.5 Gyr after that (rise is 
not so prominent with M acc ). The interpretation is that 
the cooling time (and hence min[f coo i/fff]) decreases as 
the core cools during the cooling leg of the cycle. This 
leads to the condensation of radially dominant cold gas 
after a cooling time (few 100 Myr) and an enhancement 
of the mass accretion rate and the jet power. Sudden 
increase in jet power overheats the core and min(t coo i/tg) 
increases after a lag of few 100 Myr; cold gas mass and 
M acc also decrease consequently. 

The bottom three panels in Figure 15 show that cold 
gas mass (radially dominant), mass accretion rate, and 
jet power are positively correlated. A slight skew toward 
a negative time lag shows that the cold gas mass and the 
mass accretion rate increases first, and that gives rise 
to an increase in jet power. Thus, the cross-covariance 
behavior of different variables is similar to that seen in 
cold gas-jet cycles in Figure 7. Also note in Figure 15 
that there are smaller number of oscillations for higher 
feedback efficiency (or smaller halo mass); this is a re¬ 
flection of smaller number of cooling/feedback cycles in 
these cases (see Fig. 12). 

In our 3-D simulations we see the build-up of a massive 
rotationally-supported torus. A part of this torus should 
cool further and lead to star formation, as argued by 
Li et al. (2015). While the cold gas (mostly in the 
torus) builds up in time and saturates after 2 Gyr, the 
jet energy shows fluctuations in time even after that (see 
the top panel of Fig. 9 and Fig. 14). Therefore, in our 
models there is no correlation between total cold gas mass 
and jet energy. However, there is a correlation between 
the radially-dominant cold gas mass and the jet energy 
(compare green-dotted and black dot-dashed lines in the 
top panel of Fig. 9). Thus, although most cold gas is 
decoupled from jet feedback, it is the subdominant in¬ 
falling cold gas which is powering AGN. This is in line 
with the observations of McNamara, Rohanizadegan, & 
Nulsen (2011), who find no correlation between the jet 
power and the available molecular gas. They, therefore, 
argue that most of the cold gas is converted into stars 
rather than being accreted by the SMBH. 

Finally, we compare our simulations with recent obser¬ 
vational studies of cold gas kinematics and star forma¬ 
tion. These have been studied in unprecedented detail 
in some elliptical galaxies and clusters, thanks mainly to 
ALMA and Herschel telescopes (e.g., McNamara et al. 
2014; Russell et al. 2014; David et al. 2014; Edge et 
al. 2010; Werner et al. 2014; Tremblay et al. 2012; 
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Fig. 15.— Cross-covariance of various quantities (min[£ coo l/£ff ]5 jet energy, mass in cold phase, M ac c) as a function of time lag to show 
temporal relationship between these various quantities. Cross covariance between two quantities as a function of time, as used here, is 

where —T < r < T is the time lag and 5a and 5b 


defined as: cov(a, 6; t) = ^ [5a(t + r)5b(t)dt] / 


V/o T /(T \Sb(t)\ 2 dt 


are mean-subtracted quantities. Since there is a large variation in various quantities (see Figs. 13 & 14), we take log before evaluating 
cross-covariance. For the 3-D cluster run we have used the radially-dominant cold gas mass; the cross-covariance is much weaker if we use 
total cold gas mass. 


Rawle et al. 2012). In this paper we have mainly fo¬ 
cussed time-averaged kinematics, as shown in Figures 4 
and 5. We can clearly see three kinematically distinct 
components of cold gas: a rotationally-supported mas¬ 
sive torus, ballistically infalling cold gas, and jet-uplifted 
fast cold gas. 

Observations of different clusters are snapshots at a 
particular instant, at which a particular component (e.g., 
the rotating torus, a fast outflow, or a radially distributed 
inflow; see Fig. 6) of the cold gas distribution may be 
more prominent. We will present the details of cold gas 
kinematics in various states of the ICM (with cold in¬ 
flows, outflows, and the rotating torus) in a future work. 

Some of the salient properties of the cold gas distri¬ 
bution in the fiducial run are: the rotating cold gas 
torus, when present, is more massive compared to in¬ 
falling cold gas (this component may be exaggerated in 
our simulations as we do not include star formation that 
would quickly consume some of the cold gas); the ro¬ 
tating disk rotates at the almost constant local circular 
velocity (100-200 km s” 1 for our fiducial cluster run; the 
actual value may be larger because we have ignored the 
gravitational potential due to the BCG) in form of a mas¬ 
sive torus within 5 kpc; the radially-dominant cold gas 
is much more spatially extended (out to few 10s of kpc) 


compared to the rotating torus, and the majority of this 
component also has a velocity close to the circular veloc¬ 
ity; some (about 10 5 M 0 ) radially dominant outflowing 
gas has a radial velocity as high as 1000 km s -1 (a fast 
component is seen in the observations of Russell et al. 
2014 and McNamara et al. 2014); the in-falling cold gas 
(on average) is about twice as much as the outflowing 
component. 

While the accretion rate through the inner radius 
(dominated by cold gas) is smaller than 100 Mq yr” 1 
(the accretion rate on to the SMBH is < 1% of this) at 
all times (see Fig. 9), the cooling/accretion and outflow 
rates in the cold gas can be much larger instantaneously 
because of the massive cold torus buffer (Fig. 6). 

The observations show varying cold gas kinematics 
in different systems: radially in-falling cold molecular 
clouds of 3 x 10 5 to 10 7 Mq in a galaxy group NGC 5044 
(David et al. 2014); 5 x 10 10 M 0 molecular gas pre¬ 
dominantly in a rotating disk, and about 10 1Q M 0 in a 
fast (line of sight velocity up to 500 km s -1 ) outflow in 
Abell 1835 (McNamara et al. 2014); 10 10 M 0 of molec¬ 
ular gas roughly equally divided between as rotating disk 
(velocity ~ 250 km s -1 ) and a faster (570 km s _1 ) in¬ 
falling/outflowing component in Abell 1664 (Russell et 
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al. 2014). In our simulations we observe similar compo¬ 
nents of the cold gas distribution, as shown in Figures 4 
and 5. 


5. CONCLUSIONS 

Cold-mode feedback, due to condensation of cold gas 
from the hot ICM when the local density is higher than 
a critical value (see below), has emerged as an attractive 
paradigm to interpret observations in cluster cool cores. 
In this paper we have carried out simulations of clusters 
of halo masses 7 x 10 14 M 0 and 1.8 x 10 15 M 0 with 
feedback driven AGN jets, varying the feedback efficiency 
over a large range (10 -5 — 0.01). AGN feedback is able to 
suppress cooling flows within the observational limits (by 
a factor of ~ 10) even for a feedback efficiency as low as 
6 x 10 -5 . This is the major difference from previous jet 
simulations, which use a much larger feedback efficiency 
(> 10” 3 ; Gaspari et al. 2012; Li & Bryan 2014b; Li et 
al. 2015). Because of the high feedback efficiency, the 
previous simulations attain thermal equilibrium in a hot, 
low-density core (with t coo i/tg> 10) which does not show 
cold gas and jet cycles at late times. In contrast, our low 
efficiency simulations show cooling/jet cycles even at late 
times. 

The core undergoes cooling and feedback heating cy¬ 
cles because of cold gas precipitation and enhanced ac¬ 
cretion on to the SMBH. There are more cycles for a 
lower efficiency and a larger halo mass. The cool-core 
appearance is preserved even during strong jet events. 
Even with large efficiencies, jet feedback raises the core 
entropy to several tens of keV cm 2 , and therefore cannot 
explain the non-cool-core clusters with large cores and 
entropies greater than 100 keV cm 2 . The origin of these 
non-cool-core clusters is still poorly understood (see, e.g., 
Poole et al. 2008). 

In this paper we highlight some results that were not 
emphasized in previous simulations of AGN jet feedback 
in clusters; in particular, we compare our results with 
several recent observations. Following are our major con¬ 
clusions: 

• First and most importantly, the results from differ¬ 
ent codes, different setups, and different implemen¬ 
tation of jet feedback (as long as condensation and 
accretion of cold gas is accounted for; e.g., Gaspari 
et al. 2012; Li et al. 2015) give qualitatively 
similar results. This indicates the robustness of 
the cold feedback mechanism, and the importance 
of precipitation (which occurs when t coo \/tg< 10) 
and associated feedback in regulating cluster cores. 

• We find that a feedback efficiency (defined as the 
ratio of jet mechanical luminosity and the rest mass 
accretion rate [M acc c 2 ] at ~ 1 kpc; see Eq. 6) of 
as small as 6 x 10 -5 is sufficient to suppress the 
cooling/star-formation rate in cluster cores by a 
factor of about 10 (see Fig. 8). An even smaller 
efficiency is sufficient for lower mass halos because 
the thermal energy of the ICM is smaller compared 
to the rest mass energy. Our fiducial efficiency is at 
least 20 times lower than the models of Li & Bryan 

(2014b) and Gaspari et al. (2012). Our values are 
consistent with the expectation that the mass ac¬ 
cretion rate on to the SMBH is much smaller than 


the accretion rate estimated at ~ 1 kpc, and the 
fact that powerful jets exist only when the SMBH 
accretion rate is smaller than 0.01 time the Edding¬ 
ton rate. 

The required efficiency can be roughly estimated 
as follows. On average, the core luminosity is bal¬ 
anced by the energy input rate; i.e., e req M acc c 2 ~ 
Lx (ereq is the required feedback efficiency, M acc 
is the accretion rate estimated at 1 kpc, and Lx is 
the X-ray luminosity of the cooling core). If we as¬ 
sume that the mass accretion rate is a fixed factor 
/ (<C 1) of the cooling flow value then, 


c req./ 


? M c c 2 

tcool 


1.5 k B TM c 

filTlpt COO l 


(M c is the core mass) implies that 



3 x 10 


-5 


0.1 


-1 


M, 


200 


2/3 


7 x 10 14 M 0 y 


normalizing to the parameters of our fiducial clus¬ 
ter (T = 2 keV; see the bottom right panel of Fig. 
11), where c 2 = ksT / /j,m p is the sound speed of 
the core ICM. This estimate agrees with our fidu¬ 
cial efficiency e = 6 x 10 -5 . 


• We observe cycles in jet energy, radially-dominant 
cold gas mass, and mass accretion rate which 
are governed by t coo i/tg measured in the hot 
phase. If t coo \/tg< 10 cold gas precipitates, and 
leads to multiphase cooling and enhanced accre¬ 
tion on to the SMBH. Sudden rise in the accre¬ 
tion rate, for a sufficiently high feedback efficiency, 
leads to overheating of the core and an increase in 
tcooi/^ff above the threshold for cold gas conden¬ 
sation. We emphasize that thermal equilibrium in 
cluster cores only holds in a time-averaged sense. 
There are cooling/heating cycles during which the 
core slowly cools/heats up. The core spends a 
longer time in the hot state for a larger feedback ef¬ 
ficiency and a lower halo mass, leading to a smaller 
number of cooling/heating cycles (see Fig. 12). 

Several observations hint at cycles in jet power 
and cooling of the hot gas (see Fig. 7). We do 
not expect such cycles if feedback occurs via the 
smooth hot/Bondi mode as we do not have sud¬ 
den cooling/feedback events. The hysteresis be¬ 
havior observed in the core X-ray properties (Ro, 
min[t coo i/tff]) and the mass of cold gas, jet power, 
etc. leads to a large dispersion in the correlation 
between these quantities (see Figs. 13 & 14). In 
particular, the mass accretion rate (at 1 kpc) is in¬ 
dependent of the total cold gas mass, which is dom¬ 
inated by the rotating cold torus, most of which is 
consumed by star formation (rather than accretion 
on to SMBH; e.g., McNamara, Rohanizadegan, & 
Nulsen 2011). 

• We can classify the cold gas in our 3-D simulations 
into two spatially and kinematically distinct com¬ 
ponents: a centrally concentrated (within 5 kpc), 
rotationally supported (It^l » |iy|) torus (Fig. 3); 
and extended (both infalling and outgoing) cold 
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gas going out to 30 kpc (Fig. 4). The massive, 
rotationally-supported disk is decoupled from the 
feedback loop; the radially-dominant infalling cold 
gas is what closes the feedback cycle. The cold 
torus rotates at the local circular speed (200-300 
km s” 1 ). The infalling cold gas can be fast (< 400 
km s _1 ), but the uplifted cold gas from the rotating 
torus can sometimes reach speeds larger than 1000 
km s -1 as it is accelerated by the fast jet (Fig. 5). 
The mass of the radially-dominant infalling cold 
gas is about a factor of two times the outflowing 
cold gas. The massive cold torus is expected to be 
substantially depleted by star formation, which we 
do not take into account in our simulations. 

• The minimum in the ratio of the cooling time and 
the free fall time (min[f coo i/iff]) seems better than 
the core entropy (Kq) for characterizing the cool 
cores. First of all it is a dimensionless parame¬ 


ter which applies for all halo masses (Voit et al. 
2015a), and secondly it is not sensitive to strong 
cooling or heating in the very center (unlike A'o). 
The entropy and t C ooi/te panels in Figure 11 show 
that a constant i C ooiAff ‘core’, which corresponds 
to a double power law for the entropy profile (with 
a slow increase with radius in the core; as argued 
in Panagoulia, Fabian, & Sanders 2014), is a bet¬ 
ter approximation to clusters in the very cool state 
with min(i coo i/tff)< 5. 
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